Isotope effects in ice Ih: A path-integral simulation 
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Ice Ih has been studied by path-integral molecular dynamics simulations, using the effective q- 
TIP4P/F potential model for flexible water. This has allowed us to analyze finite-temperature 
quantum effects in this solid phase from 25 to 300 K at ambient pressure. Among these effects we 
find a negative thermal expansion of ice at low temperatures, which does not appear in classical 
molecular dynamics simulations. The compressibility derived from volume fluctuations gives results 
in line with experimental data. We have analyzed isotope effects in ice Ih by considering normal, 
heavy, and tritiated water. In particular, we studied the effect of changing the isotopic mass of 
hydrogen on the kinetic energy and atomic derealization in the crystal, as well as on structural 
properties such as interatomic distances and molar volume. For D2O ice Ih at 100 K we obtained 
a decrease in molar volume and intramolecular O-H distance of 0.6% and 0.4%, respectively, as 
compared to H2O ice. 

PACS numbers: 65.40.De,71.15.Pd 
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I. INTRODUCTION 

Condensed phases of water have been studied over 
the years by many experimental and theoretical tech- 
niques. Several of their properties have been well ex- 
plained, but some others lack a complete understand- 
ing, in part due to the peculiar structure of liquid and 
solid water, in which hydrogen bonds between adjacent 
molecules give rise to properties somewhat different than 
those of most known liquids and solids (the so-called wa- 
ter "anomalies")^— 

The structure of "normal" ice Ih has been known for 
many years. A feature displayed by this and some other 
ice phases is that although water molecules are disposed 
on a regular crystal lattice, their spatial orientation is 
disordered, being only controlled by the so-called ice 
rules££ Thus, ice Ih is not an ordinary crystal, in the 
sense that oxygen atoms lay on regular lattice positions, 
but the hydrogen positions are randomly chosen among 
the six possible orientations of each molecule, in such a 
way that the orientations of neighboring molecules are 
compatible. 7 

Computer simulation of water in condensed phases at- 
tracted much interest after the development of Monte 
Carlo and molecular dynamics methods for studying liq- 
uids and solids at an atomic level. In the first works, 8,9 
rigid nonpolarizable models for the water molecule were 
employed, and subsequently a considerable amount of at- 
tention was focused on the development and refinement 
of empirical potentials to describe both liquid and solid 
phases of water. Nowadays, a large variety of this kind 
of potentials are present in the literature*^— Many of 
them assume a rigid geometry for the water molecule, 
and some others include molecular flexibility either with 
harmonic or anharmonic OH stretches. Also, in some 
cases the polarizability of the water molecule has been 
explicitly introduced into the model potentials.- 5 More- 
over, in recent years, some simulations of water using 



ab initio density functional theory (DFT) have appeared 
in the literature ! 16 ' 17 Nevertheless, the hydrogen bonds 
present in condensed phases of water seem to be difficult 
to describe with presently available energy functionals, 
which causes that some properties are poorly reproduced 
by DFT simulations. This is the case of the melting tem- 
perature of ice Ih, which can be overestimated by more 
than 100 KJ£ 

A limitation of ab-initio electronic-structure calcula- 
tions in condensed matter is that they usually treat 
atomic nuclei as classical particles, and typical quan- 
tum effects like zero-point motion are not directly ac- 
cessible. These effects can be included by using har- 
monic or quasiharmonic approximations, but are diffi- 
cult to take into account when large anharmonicities 
are present, as can happen for light atoms like hydro- 
gen. To consider the quantum character of atomic nu- 
clei, the path-integral molecular dynamics (or Monte 
Carlo) approach has proved to be very useful. A remark- 
able advantage of this procedure is that all nuclear de- 
grees of freedom can be quantized in an efficient manner, 
thus including both quantum and thermal fluctuations 
in many-body systems at finite temperatures. In this 
way, Monte Carlo or molecular dynamics sampling ap- 
plied to evaluate finite-temperature path integrals allows 
one to carry out quantitative and nonperturbative stud- 
ies of highly-anharmonic effects in condensed matterj 19 i 20 
Earlier studies of ice using path-integral simulations have 
been carried out by using mainly effective potentials, and 
were focused on structural and dynamic properties of the 
solid phase i 14 ' 21 " — 

A typical quantum effect is the isotopic dependence 
of several properties of a crystal, which would not vary 
with the atomic masses in a classical approach. Thus, 
the actual lattice parameters of two chemically identi- 
cal crystals with different isotopic composition are not 
equal, as a consequence of the dependence of the atomic 
vibrational amplitudes on the atomic mass<2£r— Lighter 
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isotopes have larger vibrational amplitudes (as expected 
in a harmonic approximation) and larger lattice parame- 
ters (an anharmonic effect). This effect is most noticeable 
at low temperatures, as the atoms in the solid feel the 
anharmonicity of the interatomic potential due to zero- 
point motion. At higher temperatures, the vibrational 
amplitudes are larger (causing a larger volume), but the 
isotope effect on the crystal volume becomes less promi- 
nent, because those amplitudes are less mass-dependent. 
In the high-temperature (classical) limit this isotope ef- 
fect disappears. Something similar happens for the inter- 
atomic distances in the solid. Path integral simulations 
have been used earlier to study isotopic effects in solids. 
In particular, this technique has turned out to be sensi- 
tive enough to quantify the dependence of crystal volume 
on the isotopic mass of the constituent atomsi2£~— 

All these quantum effects become more relevant as the 
atomic mass decreases, and will be particularly impor- 
tant in the case of hydrogen. Then, we pose the question 
of how the mass of the lightest atom can influence the 
structural properties of a solid water phase, in particu- 
lar ice Ih. This refers to the crystal volume, but also to 
the interatomic distances in the solid. Moreover, chang- 
ing the hydrogen mass should modify the actual binding 
energy of the crystal, since the kinetic energy of lighter 
atoms will be higher. 

With this purpose, we study in this paper ice Ih by 
path-integral molecular dynamics (PIMD) simulations. 
This technique allows us to analyze in a quantitative 
way several effects associated to the quantum nature of 
atomic nuclei, an in particular the influence of isotopic 
mass on the crystal volume and the interatomic distances. 
We consider normal (H2O), heavy (D2O), and tritiated 
(T2O) water. Interatomic interactions are described by 
the flexible q-TIP4P/F model, which has been recently 
developed and was employed to carry out PIMD simu- 
lations of liquid water^ In an earlier paper— we have 
used this interatomic potential to study the melting of 
ice, by using nonequilibrium techniques which allowed us 
to obtain the coexistence temperature of solid and liq- 
uid water phases at ambient pressure. Here we employ 
equilibrium PIMD to compare results for ice, obtained 
for the different hydrogen isotopes. 

The paper is organized as follows. In Sec. II, we de- 
scribe the computational method and the models em- 
ployed in our calculations. Our results are presented in 
Sec. Ill, dealing with atomic derealization, kinetic en- 
ergy, crystal volume, interatomic distances, and bulk 
modulus of ice Ih. Sec. IV includes a summary of the 
main results. 



II. COMPUTATIONAL METHOD 

A. Path-integral molecular dynamics 

Our calculations are based on the path-integral for- 
mulation of statistical mechanics. In this formulation, 



the partition function of a quantum system is evaluated 
by a discretization of the density matrix along cyclic 
paths, consisting of a hnite number L (Trotter num- 
ber) of steps i 32 ' 33 In the implementation of this tech- 
nique to numerical simulations, such a discretization 
gives rise to the appearance of L replicas (or beads) for 
each quantum particle. These replicas can be treated 
in the calculations as classical particles, since the par- 
tition function of the original quantum system is iso- 
morph to that of a classical one, obtained by replac- 
ing each quantum particle by a ring polymer made of 
L particlesi 19 ' 20 This path-integral approach allows one 
to study finite-temperature properties of quantum many- 
body problems in a nonperturbative scheme, even in the 
presence of large anharmonicitiesi 2 ^ The configuration 
space can be adequately sampled by molecular dynam- 
ics or Monte Carlo techniques. Here, we have used the 
PIMD method, mainly because in this case the codes are 
more easily parallelizable, a decisive factor for efficient 
use of modern computer architectures. 

Simulations of ice Ih have been carried out here in the 
isothermal-isobaric NPT ensemble (N, number of par- 
ticles; P, pressure; T, temperature), which allows us to 
calculate the equilibrium volume of the solid at given 
pressure and temperature. We have employed effective 
algorithms for performing PIMD simulations in this sta- 
tistical ensemble, as those described in the literaturei 34 i 35 
Sampling of the configuration space has been carried out 
at ambient pressure (P = 1 bar) and temperatures be- 
tween 25 and 300 K. To study isotope effects on the vol- 
ume and other properties of ice, we have considered H2O, 
D2O, and T2O. This is easily done in PIMD simulations, 
since the atomic mass is an input parameter in the calcu- 
lations. For comparison, some simulations of ice Ih were 
also performed in the classical limit, which is obtained in 
our context by setting the Trotter number L = 1. Also 
for the sake of comparison with the results for actual ice 
Ih, we carried out some PIMD simulations in which the 
hydrogen mass was assumed to be very large, so that for 
our purposes H behaves as a classical particle. In fact, 
we took a mass mjj = 10,000 u, which gives results indis- 
tinguishable from those obtained for even larger masses. 

Our simulations were carried out on ice Ih supercells 
with periodic boundary conditions. To check the in- 
fluence of finite-size effects on the results, we consid- 
ered two kinds of orthorhombic supercells. The smaller 
one included 96 water molecules and had parameters 
(3a, 2-\/3a, 2c), where a and c are the standard hexag- 
onal lattice parameters of ice Ih. The larger supercell in- 
cluded 288 molecules and had parameters (4a, 3-\/3a, 3c). 
Results obtained for both types of supercells coincided 
within error bars caused by the statistical uncertainty as- 
sociated to the simulation procedure. Prior to the PIMD 
simulations, proton-disordered ice structures were gen- 
erated by a Monte Carlo procedure, in such a way that 
each oxygen atom has two chemically bonded and two H- 
bonded hydrogen atoms, and with a cell dipole moment 
close to zero^ 
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For the interatomic interactions we have employed the 
point charge, flexible q-TIP4P/F model, 30 which has 
been previously used to study liquid water, the coexis- 
tence between the liquid and ice Ih phases^ as well as 
water clusters^ The Lennard- Jones-type interaction ap- 
pearing in the intermolecular potential was truncated at 
r c = 8.5 A. To compensate for this truncation, standard 
long-range corrections were computed assuming that the 
pair correlation function g(r) is unity for distance r > r c , 
leading to well-known corrections for the pressure and 
internal energy^ The electrostatic energy and associ- 
ated long-range forces were calculated by the Ewald 
method, using a parallel code to improve the speed of 
the procedure^ 

For a given temperature, a typical simulation run con- 
sisted of 4 x 10 4 PIMD steps for system equilibration, 
followed by 6 x 10 5 steps for the calculation of ensemble 
average properties. To keep roughly a constant preci- 
sion in the PIMD results at different temperatures, the 
Trotter number was scaled with the inverse temperature 
(L oc 1/T), so that LT = 6000 K, which translates into 
L — 20 and 240 for T — 300 K and 25 K, respectively. 

The simulations were performed by using a staging 
transformation for the bead coordinates. The constant- 
temperature ensemble was generated by coupling chains 
of four Nose- Hoover thermostats to each staging variable. 
To generate the NPT ensemble, an additional chain of 
four barostats was coupled to the volume^ To inte- 
grate the equations of motion, we employed a reversible 
reference-system propagator algorithm (RESPA), which 
allows one to define different time steps for the integra- 
tion of fast and slow degrees of freedom^ The time step 
At associated to the calculation of q-TIP4P/F forces 
was taken in the range between 0.1 and 0.3 fs, which 
was found to be appropriate for the interactions, atomic 
masses, and temperatures considered here, and provided 
adequate convergence for the studied variables. For the 
evolution of the fast dynamical variables, associated to 
the thermostats and harmonic bead interactions, we used 
a smaller time step St = At/4. Note that for ice at 25 K, 
a simulation run consisting of 6 x 10 6 PIMD steps requires 
calculation of energy and forces with the q-TIP4P/F 
potential for 1.4 x 10 9 atomic (classical) configurations, 
which was carried out by parallel computing. 

B. Spatial derealization 

We define here some spatial properties of the nuclear 
quantum paths that are helpful in the analysis of the 
PIMD simulation results. The center-of-gravity (cen- 
troid) of the quantum paths of a given particle is cal- 
culated as 

r = ^I>, (1) 

rj being the position of bead i in the associated ring 
polymer. 



The mean-square displacement of a quantum particle 
along a PIMD simulation run is then given by: 

A? = ±^J> - <r)) 2 ^ , (2) 

where (...) indicates a thermal average at temperature T. 
After a straightforward transformation, one can write A 2 
as 

A, 2 = Ql + Cl , (3) 

with 
and 

Cl = <r 2 ) - (r) 2 . (5) 

The first term in Eq. ©, Q 2 , is the mean-square "radius- 
of-gyration" of the ring polymers associated to the quan- 
tum particle under consideration, and gives the average 
spatial extension of the paths.— The second term on the 
r.h.s. of Eq. (J3J), C 2 , is the mean-square displacement 
of the path centroid. In the high-temperature (classi- 
cal) limit, each path collapses onto a single point and the 
radius-of-gyration vanishes, i.e., Q 2 — > 0. In cases where 
the anharmonicity is not very large, the spatial distri- 
bution of the centroid r is similar to that of a classical 
particle moving in the same potential. 



III. RESULTS 
A. Hydrogen derealization 

We first consider the spatial derealization of hydrogen 
in normal ice Ih, which is expected to include a nonneg- 
ligible quantum contribution. We have calculated sep- 
arately both terms giving the atomic derealization in 
Eq. (|3J), for each hydrogen atom in the water molecules. 
In Fig. 1 we display the values of Q 2 (spreading of the 
quantum paths, circles) and C 2 (centroid derealization, 
triangles) , as derived from our PIMD simulations for the 
H2O molecule at several temperatures. The total spatial 
derealization A 2 is shown as squares. In this plot, one 
observes that at low temperatures the quantum contribu- 
tion (spreading of the paths), Q 2 , dominates the spatial 
derealization of hydrogen, since the centroid displace- 
ment, C 2 , converges to zero as T — > K. The oppo- 
site happens at temperatures higher than 70 K, and in 
the high-temperature limit (unreachable here for stability 
reasons), the quantum contribution Q 2 would eventually 
disappear, as corresponds to the classical limit. At T 
= 250 K, close to the melting temperature predicted by 
this potential model |2l we found Q 2 = 2.85 x 10 -2 A 2 
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FIG. 1: Spatial delocalization of hydrogen nuclei (protons) in 
H2O ice Ih, as derived from PIMD simulations at several tem- 
peratures. Symbols indicate the mean-square displacement 
for the centroid C;? (triangles), and the radius-of-gyration of 
the quantum paths Q% (circles). The total delocalization 
is shown by open squares. Lines are guides to the eye. 
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FIG. 2: Mean square radius-of-gyration of the quantum 
paths, Qr, for hydrogen isotopes in ice Ih, as derived from 
PIMD simulations at several temperatures for H2O (squares), 
D2O (circles), and T2O (triangles). Error bars are in the or- 
der of the symbol size. The solid line represents for a free 
hydrogen atom. 



vs C 2 = 0.170 A 2 , i.e., the former contributes a 14% to 
the total spacial delocalization of hydrogen, A 2 . We note 
that at low temperatures the quantum paths associated 
to hydrogen have an average extension of about 0.15 A, 
much smaller than the H-H distance in a water molecule, 
thus justifying the neglect of quantum exchange between 
protons in the PIMD simulations. 

The mean-square displacement A 2 of hydrogen atoms 
derived from our PIMD simulations with the q-TIP4P/F 
potential are similar to those derived by Tanaka and 
MohantjM from classical molecular dynamics simula- 
tions with the TIP4P potential, in the temperature range 
considered by these authors (from 150 to 250 K). At 
lower temperatures, classical simulations should give val- 
ues smaller than the PIMD simulations, due to its neglect 
of atomic zero-point motion. The atomic mean-square 
displacements are usually supposed to be related with 
the melting of a solid. Thus, the Lindemann criterion 
gives a threshold for the maximum amplitude of atomic 
vibrations that can be sustained by a crystal. At the 
melting point of ice Ih we find from our PIMD simu- 
lations that this amplitude is about 6% of the H-bond 
distance, as discussed elsewhere, 31 

Quantum effects are more relevant as the mass of the 
particles under consideration becomes smaller. Thus, one 
expects that the quantum delocalization of the hydro- 
gen isotopes in ice will decrease for rising isotopic mass. 
This is shown in Fig. 2, where we have plotted the mean- 
square displacement of the quantum paths, Q 2 ., for hy- 
drogen (squares), deuterium (circles), and tritium (tri- 
angles) at several temperatures. In the zero-temperature 



limit and in a harmonic approximation, Q 2 , is known to 
scale as m^ 1 ^ 2 . At 25 K, we find a ratio Ql(R) / Ql(D) 
= 1.39, slightly smaller than the low-temperature har- 
monic expectancy of \/2. In the high-temperature limit 
Ql goes to zero, but the ratio Q 2 (H)/Q 2 (D) converges to 
the inverse mass ratio j 19 ' 41 in this case mo/n^H = 2. At 
300 K we found Q 2 (H)/Q 2 (D) = 1.73, between the high 
and low-temperature limits. For T2O we found at 300 
K, Q 2 r = 1.05 x 10~ 2 A 2 , so that Q 2 (H)/Q 2 (T) = 2.44, 
also between a ratio of y3 expected at low temperature 
in a harmonic approach, and the high-temperature limit 
given by iht/™h = 3. 

For comparison with the results of our PIMD simu- 
lations for different hydrogen isotopes in ice Ih, we also 
present in Fig. 2 the mean-square radius of gyration cor- 
responding to a free hydrogen atom (without any exter- 
nal potential, solid line). This can be calculated ana- 
lytically, as shown elsewhere ; 19 ' 41 and gives values of Q 2 . 
clearly higher than for hydrogen in ice, in the whole tem- 
perature range considered here, confirming the impor- 
tance of the hydrogen environment in its quantum delo- 
calization. This is even clearly appreciable at the melting 
temperature. 

For the centroid delocalization, C 2 , of deuterium and 
tritium in ice Ih we found values almost indistinguishable 
from those obtained for hydrogen (shown in Fig. 1) in the 
whole temperature region under consideration, which in 
turn were practically the same as the mean-square dis- 
placement of hydrogen found in classical molecular dy- 
namics simulations. 
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FIG. 3: Temperature dependence of the kinetic energy of ice 
Ih, as derived from PIMD simulations. Open symbols indicate 
results derived from simulations with full quantum motion of 
the atoms in the cell: squares for H2O, circles for D2O, and 
triangles for T2O. For comparison we also present results of 
classical molecular dynamics simulations (filled squares) . The 
filled circle corresponds to a PIMD simulation of ice with the 
hydrogen mass going to infinity. Error bars are less than the 
symbol size. Lines are guides to the eye. 



FIG. 4: Molar volume of ice Ih as a function of tempera- 
ture, as derived from PIMD simulations for H2O (squares), 
D2O (circles), and T2O (triangles). For comparison, results 
of classical molecular dynamics simulations (filled squares) 
are also shown. A filled circle corresponds to the ice volume 
obtained from a PIMD simulation with the hydrogen mass 
going to infinity. Error bars are in the order of the symbol 
size. Lines are guides to the eye. A solid diamond represents 
the experimental determination from Refs. [jjlfirjl at 273 K. 



B. Kinetic energy 

A typical quantum effect related to the atomic motion 
in solids is that the kinetic energy Ek at low temperature 
converges to a finite value associated to zero-point mo- 
tion, contrary to the classical result where Ek vanishes 
at K. Path integral simulations allow one to obtain the 
kinetic energy of the quantum particles under considera- 
tion, which is basically related to the spread of the quan- 
tum paths. In fact, for a particle at given temperature 
and isotopic mass, the larger the mean-square radius-of- 
gyration of the paths, Q 2 ., the smaller the kinetic energy. 
Here we have calculated E% by using the so-called virial 
estimator, which has an associated statistical uncertainty 
lower than the potential energy of the system i 39 ' 42 

In Fig. 3 we display the kinetic energy as a function of 
temperature for ice Ih. Open symbols indicate results de- 
rived from PIMD simulations: H2O (squares), D2O (cir- 
cles), and T2O (triangles). In each case, Ek increases as 
temperature rises, and at low temperature it converges to 
the value corresponding to zero-point motion. At a given 
temperature, Ek decreases for increasing isotopic mass. 
Thus, at T = 25 K we obtained E k = 7.59, 5.81, and 
4.98 kcal/mol for H 2 0, D 2 0, and T 2 0, respectively. The 
change in kinetic energy caused by isotopic substitution 
decreases as temperature is raised. However, we stress 
that at the melting temperature, the kinetic energy of 
deuterated ice is still more than 1.5 kcal/mol lower than 



that of normal ice, reflecting the importance of quantum 
effects in ice Ih at this temperature. 

For comparison with the results of our PIMD simula- 
tions, we also present in Fig. 3 the kinetic energy corre- 
sponding to a classical model with 37V degrees of freedom 
(N atoms): Ef — 3A^fcsT/2, solid squares. For rising 
temperature, the classical kinetic energy approaches the 
results of the PIMD simulations, in particular those cor- 
responding to the heaviest molecule, T 2 0. However, at 
300 K the classical value is still lower than £Jfc(T 2 0) by 
3.14 kcal/mol, which represents about 50% of the quan- 
tum value. To assess the importance of quantum effects 
associated to the oxygen motion, we also show in Fig. 3 
the kinetic energy obtained in PIMD simulations at 150 
K for the limit bih — > 00. In this case we obtain Ek — 
2.49 kcal/mol, i.e. 1.15 kcal/mol above the classical limit 
at this temperature. This represents about 30% of the 
increase in kinetic energy from the classical limit to T 2 
quantum ice. 



C. Crystal volume 

As commented in the Introduction, the crystal volume 
is expected to change with the isotopic mass of the con- 
stituent atoms. We have calculated the equilibrium vol- 
ume of ice Ih in our isothermal-isobaric PIMD simula- 
tions for different hydrogen masses. In Fig. 4 we present 
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the temperature dependence of the molar volume for the 
different isotopes. Results of our PIMD simulations at P 
= 1 bar are shown as open symbols: H 2 (squares), D 2 
(circles), and T 2 (triangles). One first notices that the 
change of volume with temperature is not monotonic, 
but at low temperatures it decreases as temperature is 
raised, to reach a minimum at T 100 K. At higher 
temperature, the molar volume increases, and the solid 
recovers the usual thermal expansion. Also, the volume 
of the solid decreases as the isotope mass of hydrogen 
is raised, as a consequence of the smaller spatial (quan- 
tum) derealization associated to a larger nuclear mass. 
We observe for all three isotopes the same behavior of 
the volume, with a minimum around 100 K. The volume 
changes derived from the PIMD simulations are isotropic, 
i.e., the ratio between parameters of the simulation cell is 
constant when changing the temperature, within the sta- 
tistical error bar, which agrees with experimental data. 43 
For thermodynamic reasons, one expects that in the low- 
temperature limit (T — > K) the thermal expansion of 
the solid goes to zero, i.e., dV/dT — > 0. We cannot, how- 
ever, check this point at present since PIMD simulations 
at temperatures lower than 25 K would require enormous 
computational resources. 

The molar volume (or density) of ice has been deter- 
mined experimentally in very different ways, e.g., calori- 
metric, acoustical, mechanical, x-ray, optical, or nuclear 
methods (see |3 and references therein). Measurements 
of different authors typically deviate from each other by 
up to about 0.3%, which could be caused by the den- 
sity lowering effect of aging on ice crystals or by their 
very slow relaxation to equilibrium. 44 The molar vol- 
umes 19.635 cm 3 /mol of Ginnings and Corruccini 4 ^ and 
19.633 cm 3 /mol of Dantl and Gregora 4 ^ are considered 
m Ref. the most accurate determinations at normal 
pressure and 273 K. They are plotted in Fig. 4 as a solid 
diamond. At lower temperatures, there appears some dis- 
persion in the experimental points presented by different 
authors. A minimum in the molar volume was found 
in several work a 43 ' 44 i 47 in the range between 60 and 90 
K, to be compared with the minimum at around 100 K 
derived from our PIMD simulations. The minimum mo- 
lar volume was found to be around 19.3 cm 3 /mol, some- 
what smaller than our calculated value at 100 K (19.40 
cm 3 /mol). This, along with the rather good agreement 
between calculated and experimental values at the melt- 
ing temperature, indicates that the employed interatomic 
potential underestimates the thermal expansion of the 
solid. In this context, the most remarkable point is that 
the q-TIP4P/F potential is able to capture the negative 
thermal expansion of ice at low temperatures. This phe- 
nomenon is associated to the tetrahedral coordination 
of the water molecules in ice, and has been observed in 
other solids with similar structures. In particular, it is 
well known for crystals with diamond and zinc blende 
structure, such as Si, GaAs, and CuCl4£ 

Another important point to be stressed here is the ap- 
preciable isotopic effect on the molar volume. In fact, 



at T = 100 K, we found v = 19.400, 19.285, and 19.240 
cm 3 /mol, for H2O, D2O, and T2O, respectively, which 
translates into a volume reduction of 0.6% and 0.8% 
in passing from normal ice to deuterated and tritiated 
ice, respectively. From the crystal structure obtained by 
Rottger et a/.— one derives for normal ice Ih at 100 K a 
molar volume of 19.301 cm 3 /mol, somewhat smaller than 
that found in the PIMD simulations with the q-TIP4P/F 
potential. These authors also measured the lattice pa- 
rameters of deuterated ice, and found an inverse isotopic 
effect, which is not reproduced in our simulations. 

For comparison, we also present in Fig. 4 the molar 
volume of ice Ih derived from classical molecular dynam- 
ics simulations with the same interatomic potential (filled 
squares). At a given temperature, the volume obtained 
from the classical simulations is smaller than that found 
in quantum PIMD simulations, but the difference de- 
creases as T rises. Note that, contrary to the results 
of the quantum simulations, the molar volume obtained 
in the classical ones increases monotonically and does not 
show any anomaly in the whole temperature region up to 
300 K. At T = 100 K, we obtain in the classical approxi- 
mation a molar volume v = 18.928 cm 3 /mol. This means 
that at this temperature, quantum effects cause for nor- 
mal ice Ih a remarkable volume expansion of about 2.5%. 
In the zero-temperature limit one can estimate for H2O 
ice a volume increase of about 0.8 cm 3 /mol, which means 
a rise around 4.3% due to nuclear quantum effects. 

We note that both the thermal expansion (positive or 
negative) and the isotopic effect on the crystal volume are 
anharmonic effects that would be absent in a harmonic 
approximation. On one side, the thermal expansion ap- 
pears in both classical and quantum simulations, since 
in both cases the atomic vibrations feel the anharmonic- 
ity of the interatomic potential at any finite temperature 
(T > K). In the quantum case, the anharmonicity is 
also felt in the low-temperature limit T — > K, due to 
zero-point motion, and thus in this limit the crystal vol- 
ume is larger for quantum than for classical calculations 
(which may be called "zero-point crystal expansion"). 
On the other side, the isotopic effect on the crystal vol- 
ume is a typical quantum effect that disappears in the 
classical limit, and is associated to the different vibra- 
tional amplitudes of different isotopes (which coincide in 
the classical limit). 

The negative thermal expansion of ice at low T is due 
to low-energy transverse vibrational modes with negative 
Griineisen parameter. For a mode in the n'th phonon 
branch with wave vector q, this parameter 7n(q) is de- 
fined from the logarithmic derivative of its frequency with 
respect to the crystal volume4^ 



At relatively low temperature, low-energy modes are 
more populated than modes with higher energy (with 
positive 7n(l))j and then the overall contribution to 
the volume change with increasing temperature will be 



negative! 50 ' 51 In this line, Tanaka calculated the temper- 
ature dependence of the crystal volume in a quasihar- 
monic approximation and obtained a negative thermal 
expansion at low temperatures! 50 ' 52 This author calcu- 
lated the relative weight of different vibrational modes to 
the overall thermal expansion at different temperatures, 
and found that a family of low energy modes (u> < 50 
cm -1 ), corresponding to bending motion of three water 
molecules is largely responsible for the negative thermal 
expansion at low temperatures. In connection with this, 
an instability in a transverse acoustical mode in ice Ih 
was found at low temperature from incoherent inelastic 
neutron scattering. 51 

We remember that in the classical molecular dynam- 
ics simulations of ice Ih we did not obtain any negative 
thermal expansion, in contrast with PIMD simulations. 
This is in line with the explanation given above related 
to the different weight of phonons with different energies 
at low temperatures. In a classical model all phonons 
have the same weight at any temperature (equipartition 
principle) , and thus the negative Griineisen parameter of 
low-energy transverse modes is largely compensated for 
at any temperature by the positive 7 of most vibrational 
modes in the solid. 

It is worthwhile commenting on the effect of pressure 
upon the molar volume of ice. It is well known that ice Ih 
is no longer mechanically stable at pressures in the order 
of 1 GPa, where it has been observed to transform into 
an amorphous phase i 53 ' 54 This has been in fact obtained 
in our PIMD simulations at T = 251 K and P = 1 GPa, 
where the ice crystal collapsed into an amorphous solid, 
with a volume decrease of about 18%. At P = 0.8 GPa 
and the same temperature we found that the ice crystal 
was still stable along our simulation runs, and obtained 
a reduction in the crystal volume of 8.1% with respect to 
ambient pressure. Interestingly, at 0.8 GPa the atomic 
derealization of hydrogen was found to be A 2 = 0.309 
A 2 , to be compared with a value of 0.199 A 2 obtained 
at P = 1 atm (see above), which means a large increase 
in the mean-square displacement of about 50%. Such an 
anomalous increase in the atomic derealization seems to 
be related with the crystal instability associated to the 
amorphization of ice Ih under pressure. This question 
requires further investigation by using path-integral sim- 
ulations at different pressures. 



D. Interatomic distances 

In this subsection we present results for interatomic 
distances in ice Ih between atoms in the same and adja- 
cent molecules. We first show in Fig. 5 the mean distance 
between oxygen atoms in neighboring molecules. Open 
symbols represent results of PIMD simulations: squares 
for H2O and triangles for T2O. Data points for the case 
of D2O lie between those of H2O and T2O, and are not 
shown for clarity. For comparison we also present in this 
figure the average 0-0 distance derived from classical 
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FIG. 5: Mean distance between oxygen atoms in nearest- 
neighbor water molecules, as a function of temperature. Open 
symbols correspond to results of PIMD simulations for H2O 
(squares) and T2O (triangles). Filled squares represent re- 
sults of classical molecular dynamics simulations. Lines are 
guides to the eye. 



molecular dynamics simulations (filled squares). Com- 
paring results of PIMD and classical simulations, we find 
that in the low-temperature limit the 0-0 distance for 
normal ice increases by about 0.04 A with respect to 
the classical value, i.e., around 1.5%. The 0-0 distance 
derived from quantum simulations increases as tempera- 
ture is raised, and the slope rises from a vanishing value 
at T — > K to nearly the classical value at 300 K. For 
tritiated ice we found an interatomic distance between 
that for H2O ice and the classical limit. 

For a totally homogeneous and isotropic crystal ex- 
pansion one expects a volume change AV/V = 3AZ/7, 
I being any distance in the solid. This relation is valid 
for small volume changes, i.e., for AV/V <C 1, as occurs 
in the temperature range considered here, and in par- 
ticular, one would expect it to be fulfilled for the inter- 
molecular distance in the solid. Thus, from the increase 
in the average 0-0 distance at low temperature for ice 
Ih one expects a volume expansion of about 4.5% due to 
zero-point motion. This value is close to that estimated 
from the volumes obtained in classical and quantum sim- 
ulations at low temperature, which give an expansion of 
4.3% (see above, Sec. III.C). Note, however, that a fixed 
relation between interatomic distances and crystal vol- 
ume can not be strictly followed in the case of ice Ih when 
the temperature is raised, specially in the region where 
the crystal shrinks. In fact, we observe that the average 
0-0 distance does not present a minimum atTss 100 K, 
as the molar volume of ice. This clearly reflects the fact 
that shrinking of the ice crystal is not due to a reduc- 
tion in the distance between nearest molecules, but to a 
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FIG. 6: Temperature dependence of the mean intramolec- 
ular O-H distance, as derived from simulations of ice Ih. 
Open symbols represent results of PIMD simulations for H2O 
(squares), D2O (circles), and T2O (triangles). For compari- 
son, the average O-H distance obtained from classical molec- 
ular dynamics simulations (filled squares) are also presented. 
The filled circle corresponds to the limit ran — > 00 in PIMD 
simulations. Error bars are in the order of the symbol size. 
Lines are guides to the eye. 



bending motion of contiguous molecules, corresponding 
to low-frequency vibrational modes, as indicated above 
in Sec. III.C. From the results of our quantum and clas- 
sical simulations at 100 K we find a volume increase of 
2.5% due to quantum nuclear effects vs a relative rise of 
0.9% in the average 0-0 distance, which would corre- 
spond to a volume change of 2.7% for an isotropic and 
homogeneous expansion. 

We now turn to analyze the interatomic distances in- 
side water molecules in the ice crystal, as derived from 
our simulations at different temperatures. In Fig. 6 we 
have plotted the mean intramolecular O-H distance as 
a function of temperature. As in previous figures, open 
symbols indicate results of PIMD simulations for H2O 
(squares), D 2 (circles), and T 2 (triangles). At a 
given temperature, the average O-H distance is larger for 
smaller hydrogen isotopic mass. In particular, the differ- 
ence between O-H and O-D distances amounts to about 
4 x 10~ 3 A in the whole temperature range under consid- 
eration, which agrees with data derived from Compton 
scattering experiments for water at room temperature. 55 
For comparison, we present in Fig. 6 results for the clas- 
sical limit (filled squares), which lie clearly lower than 
those of the quantum simulations. To assess the rela- 
tive importance of quantum nuclear effects of hydrogen 
and oxygen, we give also the O-H distance obtained in a 
PIMD simulation in the limit of infinite hydrogen mass 
(filled circle at 150 K). This point is still clearly above 
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FIG. 7: Temperature dependence of the average intramolec- 
ular O-H distance, as derived from PIMD simulations of a 
single H2O molecule (triangles), ice Ih (squares), and liquid 
water (circles). Error bars are in the order of the symbol size. 
Lines are guides to the eye. 



the classical limit, with a quantum contribution (coming 
from quantum derealization of oxygen nuclei) amount- 
ing to 31% and 20% of the increase in O-H distance for 
full quantum T 2 and H 2 ice, respectively. 

It is interesting that the intramolecular O-H distance 
decreases for increasing temperature, in both classical 
and quantum approaches. This fact has been observed 
experimentally and reported in the literature.— It is due 
to the peculiar structure of ice with hydrogen bonds con- 
necting neighboring water molecules. As temperature is 
raised, molecular motion is enhanced, so that hydrogen 
bonds become softer, and the average intermolecular O- 
H distance rises (in agreement with the increase in 0-0 
distance presented above). This causes a strengthening 
of the intramolecular O-H bonds, with a concomitant 
decrease in the interatomic distance in water molecules. 
The same behavior is found from the classical molecular 
dynamics simulations (filled squares in Fig. 6), although 
in this case the mean intramolecular O-H distance is 
clearly smaller than in the full quantum ice. 

To obtain deeper insight into the contraction of the in- 
tramolecular O-H bond for rising temperature, we com- 
pare in Fig. 7 the average O-H distance derived from 
our PIMD simulations for ice Ih (open squares) with 
those obtained for liquid water (circles) and for an iso- 
lated H 2 molecule (triangles). One observes first that 
the O-H distance in the water molecules increases ap- 
preciably in passing from the gas phase to condensed 
phases (solid and liquid). This is a consequence of hy- 
drogen bonds with the nearest molecules, which weaken 
the intramolecular bonds. The enlargement of these in- 
tramolecular bonds is larger for the solid than for the 



liquid, since in the former the intermolecular hydrogen 
bonds are stronger due to the lower mobility of the 
molecules and the larger rigidity of the H-bond network. 
In the case of ice Ih, the bond enlargement with respect 
to an isolated H2O molecule amounts to about 0.03 A, 
which means a 3% of the bond distance. Second, one 
observes that the O-H distance for an isolated molecule 
increases slightly as the temperature is raised, contrary 
to the condensed phases, where this distance decreases 
for increasing temperature. The behavior observed for 
the isolated molecule is the usual thermal expansion due 
to anharmonicities in the interatomic potential. For the 
condensed phases, this natural expansion is largely com- 
pensated for by intermolecular interactions through hy- 
drogen bonds, as indicated above. As a result, the in- 
tramolecular O-H distance in these phases decreases as 
the temperature rises. 



E. Bulk modulus 

Among the many surprising properties of liquid water 
and ice, one finds that their compressibility is smaller 
than what one could naively expect from the large cav- 
ities present in their structure, which could presum- 
ably collapse under pressure without water molecules ap- 
proaching close enough to repel each other. For ice Ih, in 
particular, the hydrogen bonds holding the crystal struc- 
ture are rather stable, as manifested by the relatively 
high pressure necessary to break down the H-bond net- 
work and amorphize the solid^ 

The isothermal compressibility k of ice, or its inverse 
the bulk modulus [B = 1/k — -V{dP/dV) T ] can be 
calculated directly from our PIMD simulations in the 
isothermal-isobaric ensemble. In fact, in the NPT en- 
semble the isothermal bulk modulus can be obtained 
from the mean-square fluctuations of the volume, a v — 
(V 2 ) — (V) 2 , by using the expression 57,58 



B = 



k B T(V) 



(7) 



where ks is Boltzmann's constant. This expression 
has been employed earlier to calculate the bulk mod- 
ulus of different types of solids from path-integral 
simulations! 58 ' 59 Here, we have used Eq. (J7J to obtain 
B for ice Ih from both classical and PIMD simulations. 

In the classical simulations we find that the bulk mod- 
ulus decreases roughly linearly as temperature is raised, 
and in the zero-temperature limit it extrapolates to a 
value Bq = 17.6 ± 0.3 GPa. In the quantum simulations 
of H2O ice Ih, the solid is found to be "softer" than in 
the classical simulations, in the sense that in the former 
the volume fluctuations are larger, and consequently the 
bulk modulus becomes smaller. Note that in Eq. {7} one 
also has the average volume (V") in the numerator, which 
is larger for the quantum than for the classical solid (see 
Sec. III.C), but this volume increase due to quantum nu- 
clear effects is dominated by the also present larger vol- 
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FIG. 8: Bulk modulus of ice Ih as obtained from quantum 
PIMD (open squares) and classical molecular dynamics (open 
circles) simulations for H2O. Error bars show the statistical 
uncertainty in the values of B found from the simulations. 
The filled square corresponds to the value obtained at 273 K 
from the equation of state derived by Feistel and Wagner 44 
from experimental data. Lines are guides to the eye. 



ume fluctuations a 2 ?. In the low-temperature limit we 
find an extrapolated value of B = 13.8 ± 0.4 GPa, about 
3.8 GPa smaller than the classical value, which means 
that the quantum effect reduces the bulk modulus by 
more than 25%. 

Note that the values of B derived from our simulations 
in the isothermal-isobaric ensemble show relative error 
bars larger than those corresponding to other calculated 
variables (e.g., kinetic energy, interatomic distances, or 
molar volume). This is a consequence of the statistical 
uncertainty of the volume fluctuations cry, that are em- 
ployed to calculate the bulk modulus B. The error bars 
are similar for both classical and PIMD simulations. 

There appears in the literature a dispersion of data for 
the isothermal compressibility (or bulk modulus) of ice 
Ih at the melting temperature and normal pressure. The 
bulk modulus of ice found from our PIMD simulations 
is close to that derived by Feistel and Wagner^ who 
obtained an equation of state for ice Ih from a set of 
experimental data. This value at 273 K is shown in Fig. 8 
as a solid square. It is hard to estimate the precision of 
this value of the bulk modulus due to lack of error bars 
for the published data, and the uncertainty may be large, 
as suggested by the dispersion in the experimental data 
obtained by different authors. The uncertainty is even 
larger at lower temperatures, where we could not find 
direct experimental data. From the equation of state by 
Feistel and Wagner— one derives at low temperature and 
normal pressure an isothermal bulk modulus of 10.6 GPa, 
somewhat lower than the low-temperature extrapolation 
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of our simulation results. 

The bulk moduli for D 2 and T 2 ice derived from 
our PIMD simulations show a similar trend to that of 
normal ice Ih, and have not been presented in Fig. 8 for 
clarity. They seem to lie between the classical limit and 
the quantum value for H2O, but due to the data disper- 
sion, the difference between B values for the different hy- 
drogen isotopes can be hardly quantified. Moreover, we 
do not observe in the bulk modulus of ice any anomaly 
similar to that of the crystal volume at T w 100 K, but 
we cannot exclude it from our present results because of 
their statistical uncertainty. 

IV. SUMMARY 

We have presented results of PIMD simulations of ice 
Ih with different hydrogen isotopes in a wide range of 
temperatures. This technique allows us to explore readily 
isotope effects, since the atomic masses appear as input 
parameters in the calculations. This kind of simulations 
enable one to calculate separately kinetic and potential 
energies at finite temperatures, taking into account the 
quantization of nuclear motion. This includes considera- 
tion of zero-point motion of the atoms in the solid, which 
can be hard to analyze by analytical calculations in the 
presence of light atoms and large anharmonicities. 

This kind of quantum simulations are necessary to re- 
produce some observed facts of ice, that are not captured 
by classical simulations. In this sense, PIMD simulations 
with the q-TIP4P/F potential are able to reproduce the 
negative thermal expansion of ice Ih at low temperatures. 
Also, these simulations reproduce the apparently anoma- 
lous decrease of the intramolecular O-H distance for in- 
creasing temperature. A good check of the employed in- 



teratomic potential is the calculation of a macroscopic 
observable such as the bulk modulus, which has been ac- 
curately reproduced in our calculations. Note that the 
use of a flexible potential model for water has allowed us 
to look in a realistic way at changes in the intramolecular 
O-H distance and quantum derealization of hydrogen at 
different temperatures. 

PIMD simulations give reliable predictions for several 
isotope effects in condensed matter. For ice Ih, in par- 
ticular, we have seen how the isotopic mass of hydrogen 
affects the kinetic energy and atomic derealization in 
the crystal, as well as the interatomic distances and mo- 
lar volume. Thus, we found for D2O ice Ih at 100 K a 
decrease in the crystal volume and intramolecular O-H 
distance of 0.6% and 0.4%, respectively, as compared to 
H 2 ice. Our simulations, however, do not reproduce the 
inverse isotopic effect in the crystal volume, observed in 
Ref. 43. 

Simulations similar to those presented here can be car- 
ried out for ice under an external pressure, and in partic- 
ular to study the amorphization process of the ice crystal. 
Also, the question of quantum tunneling and diffusion of 
hydrogen in the different ice structures needs further in- 
vestigation. These are challenging problems that could 
be addressed in the near future by using path-integral 
simulations. 
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